home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Libris Britannia 4
/
science library(b).zip
/
science library(b)
/
SCIENTIF
/
0807.ZIP
/
KEPLER.BAS
< prev
next >
Wrap
BASIC Source File
|
1987-09-30
|
5KB
|
211 lines
' PROGRAM "KEPLER"
' NUMERICAL SOLUTION OF KEPLER'S EQUATION FOR
' ELLIPTIC, PARABOLIC AND HYPERBOLIC ORBITS
' PUBLIC DOMAIN
' IBM-PC << QUICKBASIC COMPILER VERSION 3.0 >>
DEFDBL A-Z
CONST PI=3.14159265D0
CONST RTD=180/PI
CONST DTR=PI/180
COMMON SHARED ECC,MANOM.RAD,EANOM.RAD
'**********************************************************
CLS
PRINT
PRINT
PRINT "Program KEPLER"
PRINT
PRINT "Microsoft QuickBASIC Compiler"
PRINT "(C) Copyright Microsoft Corp. 1982-1987"
PRINT
CALL KEYCHECK
CLS
PRINT
PRINT
PRINT "Program introduction ( y = yes, n = no )"
INPUT A$
IF INSTR("yY",A$) THEN CALL INTRODUCTION
DO
CLS
PRINT
PRINT
PRINT TAB(15); "Kepler's Equation"
PRINT
PRINT
PRINT "Mean anomaly ( degrees [ 0 - 360 ] )"
INPUT MANOM.DEG
MANOM.RAD=DTR*MANOM.DEG
PRINT
PRINT "Orbital eccentricity ( non-dimensional )"
INPUT ECC
IF ECC<>0D0 THEN
CALL KEPLER
ELSE
EANOM.RAD=MANOM.RAD
END IF
' PRINT RESULTS
FORMAT$="####.######"
CLS
PRINT
PRINT
PRINT TAB(22);"Kepler's Equation"
PRINT
PRINT
PRINT TAB(5);"Eccentricity ( non-dimensional )";TAB(45);
PRINT USING FORMAT$;ECC
PRINT
PRINT TAB(5);"Mean anomaly ( degrees )";TAB(45);
PRINT USING FORMAT$;MANOM.DEG
PRINT
IF ECC<1 THEN PRINT TAB(5);"Eccentric anomaly ( degrees )";TAB(45);
IF ECC=1 THEN PRINT TAB(5);"Parabolic anomaly ( degrees )";TAB(45);
IF ECC>1 THEN PRINT TAB(5);"Hyperbolic anomaly ( degrees )";TAB(45);
PRINT USING FORMAT$;EANOM.RAD*RTD
CALL KEYCHECK
CLS
PRINT
PRINT
PRINT "Another selection ( y = yes, n = no )"
INPUT SELECTION$
LOOP UNTIL INSTR("nN",SELECTION$)
END
'**********************************************************
SUB KEPLER STATIC
' KEPLER'S EQUATION SUBROUTINE
EANOM.RAD=LOG(2*MANOM.RAD/ECC+1.8D0)
IF ECC<1D0 THEN CALL ESTIMATE
DO
IF ECC<1D0 THEN L=ECC*SIN(EANOM.RAD)
IF ECC>=1D0 THEN L=.5D0*ECC*(EXP(EANOM.RAD)-EXP(-EANOM.RAD))
IF ECC<1D0 THEN M=EANOM.RAD-L-MANOM.RAD
IF ECC>=1D0 THEN M=L-EANOM.RAD-MANOM.RAD
IF ABS(M)<=1D-8 THEN EXIT DO
IF ECC<1D0 THEN P=ECC*COS(EANOM.RAD)
IF ECC>=1D0 THEN P=.5D0*ECC*(EXP(EANOM.RAD)+EXP(-EANOM.RAD))
IF ECC<1D0 THEN Q=1D0-P
IF ECC>=1D0 THEN Q=P-1D0
R=-M/Q
R=-M/(Q+.5D0*R*L)
R=-M/(Q+.5D0*R*L+R*R*P/6D0)
R=-M/(Q+.5D0*R*L+R*R*P/6D0-R*R*R*L/24D0)
EANOM.RAD=EANOM.RAD+R
LOOP
END SUB
'**********************************************************
SUB ESTIMATE STATIC
' EMPIRICAL ESTIMATE SUBROUTINE
A=ECC^2
B=MANOM.RAD^2
EANOM.RAD=A*(-.0579437D0*B+.02665324D0*MANOM.RAD+.05868658D0)
EANOM.RAD=EANOM.RAD+ECC*(-.19174923D0*B+.46905672D0*MANOM.RAD+.600725329D0)
EANOM.RAD=EANOM.RAD-.04987627D0*B+1.17552263D0*MANOM.RAD-.12324871D0
END SUB
'**********************************************************
SUB INTRODUCTION STATIC
' INTRODUCTION SUBROUTINE
CLS
PRINT
PRINT TAB(10);"PROGRAM 'KEPLER' IS AN INTERACTIVE BASIC"
PRINT TAB(10);"PROGRAM WHICH CAN BE USED TO SOLVE"
PRINT TAB(10);"KEPLER'S EQUATION. THE ORBITAL MOTION"
PRINT TAB(10);"OF SATELLITES, COMETS AND ALL CELESTIAL"
PRINT TAB(10);"BODIES CAN BE DETERMINED FROM KEPLER'S"
PRINT TAB(10);"EQUATION. THIS EQUATION IS WRITTEN AS"
PRINT
PRINT TAB(10);" M = E - e SIN E "
PRINT
PRINT TAB(10);"FOR ELLIPTIC ORBITS OR"
PRINT
PRINT TAB(10);" M = e SINH F - F"
PRINT
PRINT TAB(10);"FOR PARABOLIC OR HYPERBOLIC ORBITS."
PRINT
PRINT TAB(10);"KEPLER'S EQUATION IS A TRANSCENDENTAL"
PRINT TAB(10);"EQUATION WHICH IS USUALLY SOLVED WITH"
PRINT TAB(10);"AN ITERATIVE NUMERICAL METHOD."
CALL KEYCHECK
CLS
PRINT
PRINT TAB(10);" IN BOTH OF THESE EQUATIONS, 'M'"
PRINT TAB(10);"IS CALLED THE 'MEAN ANOMALY' AND 'E'"
PRINT TAB(10);"IS CALLED THE 'ORBITAL ECCENTRICITY'."
PRINT TAB(10);"IN THE FIRST EQUATION, 'E' IS CALLED"
PRINT TAB(10);"THE 'ECCENTRIC ANOMALY' WHILE IN THE"
PRINT TAB(10);"SECOND EQUATION 'F' IS CALLED THE"
PRINT TAB(10);"'PARABOLIC' OR 'HYPERBOLIC' ANOMALY."
PRINT
PRINT TAB(10);"THE TRIG FUNCTION 'SINH' IN THE SECOND"
PRINT TAB(10);"EQUATION IS THE HYPERBOLIC SINE."
PRINT
PRINT TAB(10);" THE TYPE OF ORBIT IS DETERMINED"
PRINT TAB(10);"BY THE ECCENTRICITY."
PRINT
PRINT TAB(10);" 0 < e < 1 ELLIPTIC"
PRINT TAB(10);" e = 1 PARABOLIC"
PRINT TAB(10);" e > 1 HYPERBOLIC"
CALL KEYCHECK
CLS
PRINT
PRINT TAB(10);"THE ACTUAL POSITION OF A BODY IN ITS"
PRINT TAB(10);"ORBIT IS CALLED 'TRUE ANOMALY' AND"
PRINT TAB(10);"CAN BE DETERMINED FROM ECCENTRICITY"
PRINT TAB(10);"AND THE ECCENTRIC OR PARABOLIC OR"
PRINT TAB(10);"HYPERBOLIC ANOMALY."
CALL KEYCHECK
END SUB
'**********************************************************
SUB KEYCHECK STATIC
' CHECK USER RESPONSE SUBROUTINE
PRINT
PRINT
PRINT TAB(15);"< press any key to continue >"
A$=""
WHILE A$=""
A$=INKEY$
WEND
END SUB